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The principal problem considered is that of determining which placement of ?t disks 
of equal radius will cover as much as possible of a circular area A. Extensive computer 
experiments were performed to find the optimal arrangements and to compare the perform- 
ances of several ''black box" maximization methods as applied to this problem. A second 
version, in which A is divided into subregions and each disk is regarded as contributing to 
the coverage of only one subregion, is also treated. Related mathematical results and 
ciuestions are discussed. 



1. The Problem ' 
1.1. Description of Covering Problem 

We are given a circular area A of radius R, cen- 
tered at the origin of the X) '-plane, and a specified 
number n of circuhn* disks d (l<i<n) all having tlie 
same radius r<^R. How slioidd the n disks be placed 
so that they cover as great a portion of A as possible? 
And for this optunal placement of the disks, what is 
the ratio between (1) the area of the portion of A 
covered by the disks, and (2) the total area of A? 

To describe the problem more precisely we shall 
make the following defiiutions. A placement or con- 
figuration of n disks is unicjuely detennined by spec- 
ifying the (X, F) coordinates of centers of the disks. 
If we suppose that the n disks are ordered by their 
indices 6^1,62, • • • , O^, we may construct the vector 
X=(X„J\,X2,J% . . ., X„FJ where (^,,F,) de- 
notes the center oi the disk Cf. This vector X with 
2u components completely determines a configuration 
of the n disks with i-espect to the large area A. The 
configuration of figure 1 would be represented by 
the vector. 

X=(- 11,6,- 1,12,-3,3,-5,-4,6,1,12,7). 

Now for each vector there is a uniquely determined 
area of that region of the plane which is within A and 
at least one of the disks Cf. In set theoretic notation 



this region would be given as 



(,y,^'> 



A. In figure 



1 the region we want is shaded. It should be clear 
that the ratio referred to above is a function F(X). 
If we restrict tlie pairs (X,,F,) by requiring that the 
centers of all disks lie within A, tlien our problem is 
that of maximizing the function F(X) over some 
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1.2. Some Related Problems 

There are several problejus that are closely related 
to oin- problem and wdiicli seein to be somewhat more 
interesting from a purely mathematical viewpoint. 
Since their solution for the most part involves ob- 
taining a solution to the general problem stated in 
section 1.1, we shall briefly mention these others: 

(a) Given disks of radius r<^Ii, what is the mini- 
mum number of them rcquircnl to cover A com- 
pletely? 

(b) Given r<^R, what is the ]naximum number of 
disks of radius r that can be packed into A so that 
there is no overlap between disks and (^ich disk lies 
entirely within yl? See figure 2. 

(c) Given the number n of disks, wliat is the mini- 
mum radius /' for which these disks can completely 
cover Al 

(d) Given the nund)er n of disks, what is the jnaxi- 
mum I'adius r consistent with packing? 

(e) Given the value of P=nr\ what values of n 
and /' determine the best coverage? 

Certain of the above classes of problems lend them- 
selves to direct analytical solution. For exajnple, 
problem (a) above with ^ r=l/2 can be solved briefly 
as follows : For complete coverage it is required that 
the circumference of A be covered. Remembering 
that a regular hexagon inscribed in a circle of radius 
R has edges of length 7?, we see that at least six disks 
are required to cover the circumference. But if 
exactly six are used then the center of A is left un- 
covered and a seventh disk is recjuired. It can then 
be shown that seven disks are in fact suflftcient and 
the problem is solved. A shnilar argument can be 
used to show that r=l/2 is the solution to (c) when 
n=l . Problem (b) in the case r=l/2 can also be 
solved easily. 

Probably tlie least trivial analytical solution in 
this class of problems is due to Neville [1] ^, who 
solved problem (c) in the case n=5. He showed 



2 We assume /? = ! for I'onvenience; it should 1ie noted that tlie ratio of coverage 
depends only on r//? and theicfore our assumption involves 110 loss of ficiierahty. 

3 Figures in brackets indicate the Hteraturc references at the end of this paper. 
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Figure 1. Coverage by six disks. 




Figure 2. Loose packing of three disks. 



that the minimum radius required is approximately 
0.609. It is interesting to note that the configura- 
tion of disks that achieves complete coverage with 
this radius does not have central symmetry. In 
fact the boundaries of three disks (see fig. 3) pass 
thru a point near the center of A, whereas the 
other two disks are considerably displaced from its 
center. If three of the disk-boundaries are required 
to pass through the center of A, the minunum radius 
needed for complete coverage increases to 0.610, and 
it rises to 0.618 if all five boundaries are required 
to pass through the center. 

Neville^s result indicates that our intuitive expecta- 
tions, concerning the symmetr}^ of solutions of such 
problems, are not necessarily reliable. Accordingly 
no symmetry conditions were pre-imposed in the 
following work. It would be interesting to investi- 
ate further what symmetrj^ properties can be 
asserted for the configuration maximizing F(X), 
and for the configurations yielding ^ 'local maxima." 
It would be quite helpful, for possible subsequent 




Figure 3. Neville's five disk covering. 
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Figure 68. Computer configuration for Neville case. 



research, to be assured that the maximizing con- 
figurations within certain restricted classes of 
symmetric patterns actually do represent at least 
local maxima for the covering problem. Many of 
the accompanying drawings (discussed later in the 
text) display a high degree of regularity, and the 
deviations from symmetry may well be due to our 
use of a discrete grid (see sec. 3), as well as the 
inevitable inexactitude of draftmanship. 

1.3. Some Related Mathematical Literature 

The several paragraphs tliat follow contain 
references to some mathematical papers that are 
relevant to the problems discussed in this paper. 
We hope the interest of these topics will be an 
adequate compensation for the lumpiness of their 
presentation. 

An interesting result concerning coverage by disks 
of equal size is the following theorem of R. 
Kerschner [2]. 
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If N(r) is the smallest number of disks of radius 
r needed to cover a plane set of area A, then 

lim7r/'W(r)=27rV3 A/9, 

This formula suggests the possibility of using 

Aro(r)-27rV3/9r2 

as an estimate for the minimum nuinber of disks of 
radius r required to cover completely a large disk 
whose area is w (i.e., whose radius is K=l). For 
example, A^o (3/8) =9 to the nearest integer whereas 
10 disks at this radius can be made to cover 99 
percent of the area A, as indicated in table 1. 

The following theorem of Verblunsky [3] relates to 
how fast the convergence is in Kerschner's result: 

There is a number c>l/2 such that, for all small 
enough /' 

Ar(r)-(2V3/9/-^)>(2v3c/9/-), 

where N(r) is the least number of disks of radius r 
required to cover a square of area 1. This means 
that the approximation A^o(^) suggested above con- 
verges at best on the order of (1/r). In the absence 
of other knowledge, however, this might be used to 
get some idea of what sort of coverage might reason- 
ably be expected with a particular pair (n,/). The 
Verblunsky result applies to coverage of a square 
but it seejiis quite likely that the convergence is 
sunilar for the circular coverage pro])lein. 

The following result [4] establishes a relationship 
between the problems of packing and covering. If 
Ti and /'2 denote respectively the maxuimm r for 
packing and the minimum r for covering with // 
disks then 

3r2>4ri. 

This is true whenever the region to be covered or 
packed is convex. We can use this result to get a 
lower bound on 72 if we know ri, and vice vei'sa. 

The two-dhnensional case of a more general 
theorejn of D. Gale [5] implies that any plane set of 
diameter 2 can be covered by three properly chosen 
sets each of diameter <V3. The author also points 
out that no three sets each of diameter <C^J^ will 
cover a disk of diameter 2. This essentially solves 
proble^n (c) for the case n=3 and indicates that a 
disk is the '^hardest^^ set to cover among sets of 
equal diameter. 

An interesting related pi-oblcMu arises if the n disks 
d are '^thrown down'' in(h^pendentty at random, i.e., 
with their centers uniforndy distributed subject only 
to the condition that they overlap the circular area A. 
Here F(X) l)ecomes a random variables whose mean 
is an appropriate reference point in deciding which 
values of F(X) inight be considered high ones. The 
value of such a ic^ference ])oint is enhanced if one 
also has at hand the standai'd deviation a of F(X), 
which can be obtained from 



a2=M2-(M,y 

where M2 is the second movement of F(X) and Mi 
is the mean (i.e., the first moment). 

Such problems of ^ ^random coverage" have been 
treated in the technical literature. The basic 
theorem on this sul)ject, due to A. Kohnogoroff [6] 
can be stated for our purposes as follows: Suppose 
one has a probability distribution over a specified 
class of sets S in m-dijnensional Euclidean space.'* 
Then the measure ^ jLt(*S') of a set S is a random 
variable. If points of tlie Euclidean sj)ac(^ are 
denoted x=(j^i, . . ., x,J and?/=(?/i, . . ., i/j, then 
the mean of iJi(S) is given by ^ 

Mi= Prob (x is in S) dxi . . . dxm, 

its second moment by 

Al2= Prob (x is in S and y is in S) 

dxi . . . dx,r4yi . . . dyni, 

and siniil;n"l\' for highei* moments. 

This tlieorem was rediscovered by H. E. Robbins 
[7], who used it to study the one-dhnensional analog 
of our problem, i.e., random coverage of a linear 
interval by smaller intervals. He calculated Mi and 
M2 for this case, and observed that his formula for 
Ml renniins valid for the two-dimensional ('^circles") 
case that concei'ns us here. Subsetiuently J. Bronow- 
ski and J. Neynian [S] treated the randojii covei'age 
of a fixed rectangk^ by smaller rectangles with sides 
parallel to those of the fixed one. Robbins [9| 
solved the m-dimensional generalization of the 
problem for rectangles, and also treated random 
coverage of a rectangle by circular disks. L. A. 
Santalo [10] treated random coverage of an m- 
dhnensional rectangle by spheres,^ and also the 
coverage of a (two-dimensional) rectangle by rec- 
tangles of random orientation. He also solved the 
problem of random coverage of a sphere in m- 
dimensional space by sinaller splieres, which for 
m = 2 is the problem that concerns us. 

Of the many formulas derived in these papers, 
only two will be cited here. Both refer to the area 
of a circular region A of radius R=l, which is cov- 
ered by the union of n circular disks Ci of radius r<Cl, 
whose centers are independently cliosen and uniformly 
distributed over a disk of radius l+r concentric 
witli A. 

The first formula, due to Robbins, gives the mean ^ 
of this '^random covered area" as 

•''-'['-('-(.tJ)"]^ 



^ 111 our case m=2 and the sots S are non-empty intersections of the circular 
disk .1 with the union oin circular disks of radius r. 

'-■ "Measure"' is here a generic term which means "length" in one-dimensional 
situations, "area" in two dimensions, and "volume" in three. 

6 The integral formally extends over the entire Euchdean space, but in most 
applications the integrand is zero outside some bounded region. 

' Note that a "sphere" is just a linear interval in one-dimensional situations, 
and is a circular disk in two dimensions. 

* See table for values of A/i/tt pertinent to this study. 



183 



this must be divided by the area ir of A to obtain 
the mean of the ratio F{X). The second formula, 
due to Santalo, gives the corresponding variance as 



"-p^ 



-27rr^-2r^ arc cos {t/2r)+it{4r^-t''y 
7r(l+r)2 



(2 arc cos (t/2))-it(4:-t^y^^]t dt 

+^l-TY^^2)"{^'-27r(2(2r2-l) arc cos r 

-3r(l-r2)^/2_^7r+2r(l-r2)3/2 

-arcsinr)}-7r2(l-(^JJ'' 

this must be divided bv w^ to obtain the variance 

of F(X). 

2. Attempts at Analytical Solution 

2.1. Formula for Maximand 

Returning to the main problem presented in sec- 
tion 1.1, we shall describe some attempts that were 
made to obtain an analvtical solution. Problems 




Figure 4. Parameters in the two disk case. 
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Figure 5. Approximating the coverage. 



involving the maximization of a function of several 
variables can usually be handled by calculus if the 
function F(X) can be written as an expression in- 
volving the components of X and familiar functions 
of them. The first step in any attempt at an ana- 
lytical solution to our problem is to obtain some 
^ 'formula ^^ for that portion of the area of A that is 
covered by the configuration X= {Xi,Yi, . . ., Xn,Yn). 

2.2. The Two-Disk Problem 

An initial attempt was made to derive a ' 'formula'' 
for the area covered b}^ two disks of radius r<^R=l. 
The parameters describing the placement (see fig. 4) 
of the disks were 

(1) di and d2, the distances from the center of A 
to the centers of disks Ci and C2, respectively, 

(2) 6, the angle between these two distances (d<T). 
It was thought that with these parameters in place of 
{Xi,Yi,X2,Y2), it would be easier to obtain the for- 
mula desired.^ It was found that a single formula 
could not be obtained for the area covered but an 
algorithm was devised which uses no less than eight 



9 Every configuration of n dis^s can actually be specified by only (2??— 1) 
variables, by arbitrarily setting Xx=Q. This involves no loss in generality, for 
if we are given a configuration where Xij^O, then a rotation of the coordinates 
can be performed so as to make Xi =0. Such a rotation will not alter the coverage 
of the configuration. 
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Figure 6. Reflection short-cut. 
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^^forjiuilas," (Icpciuliiig on certain gcoiiietric proper- 
ties of the covering configuration. ^° 

2.3. Abandonment of Analytical Methods 

The impossibility of obtaining any reasonable 
'^formula" for the function we are trying to maxijnize 
in the relatively trivial case 71=^2 seems to indicate 
the futility of the analytical approach especially 
when n is larger. On this sad note the general 
analytical approach was abandoned and anotlier 
method of a somewhat experimental nature, using 
high-speed electronic computers, was adopted. 

3. Black Box Maximization 

3.1. General Description 

Procedures collectively knowji as '^Black Box 
Maximization' ' have been used recently to search for 
the maximum value of a function [11-14]. They are 
employed when the following conditions exist: 

(1) It is required to maximize a certain function 
F(X), where X ranges over some finite set of 
objects S. 

(2) For each individual ^it is possible to calculate 
F(X), but there is no neat analvtical expression for 
F{X}. 

(3) The function F lias some sort of continuity 
which makes it possible to define, for each X in S, a 
subset N(X) of the points of S in such a way that 
the value of F eit X differs from the value of F at 
any point of N(X) by a small amount. This set 
N(X) is called the set of neighbors of X. 

In one such method, the search for the maximum 
of F over S proceeds as follows : Pick some point Xi 
in S as a starting point and calculate F(Xi). Then 
determine the points belonging to N'(Xi) and com- 
pute F for each of these points. Among the members 
of N(Xi) select one that yields tlie highest function 
value. If this point X2 has a higher function value 
than Xi, then repeat the process with X2 replacing Xi. 
However, if F{X2)<F(Xi), then Xi is a relative 
maximum of F and the iteration terminates. 

The method we have just described is known as 
the method of ^ ^steepest ascent, '^ shice in selecting 
the new point X2 we picked that nu^mber of N(Xi) 
with the greatest function value. Two other 
methods which might be employed are worthy of 
notice. If we select as X2, that neighbor of Xi 
which has the smallest value of F among those whose 
function values are greater than F{Xi), this is known 
as the method of '^slowest ascenf' or ^^least positive 
ascent.'^ This method derives its rationale by 
analogy with the case of searching for the absolute 
maximum of a function of 2 variables where the 
function can be considered as a surface in 3 space 
witli hills and valleys representing extrema. There 
is some intuitive reason [11] for sunnising that fol- 
lowing a ^ 'river bed^^ may lead to a higher peak than a 
^^steep climb'' would obtain. The third method, 
known as the '^first positive ascent," derives its 



See appendixes 1 and 2. 



valu(* from the fact that the selection of the new 
point X2 at each stage generally takes less computa- 
tion than for the other two methods, and thus saves 
valuable tune when an electronic computer is being 
used to solve the problem. In this method the 
members of N(Xi) are arbitrarily ordered into a 
sequence Ni,N2, . . ., Np. F(Ni) is computed and 
if it exceeds F(Xi) then A^i is selected as X>. If 
F(Ni)<F{Xi) the computation is repeated with ^V2 
and so on until one of the neighbors of A"i is selected 
or all the neighbors are exhausted. W tlie latter 
occurs, then no points of N(Xi) have a higher func- 
tion value than Xi, so Xi is a relative maxijnum of F. 
There is no known way of selecting one of these 
methods as best, even given certain characteristics 
of the maximand. Gleason [11] presents interesting 
statistics comparing the ^ ^steepest ascent" with the 
^^slowest ascent" for one particular problem, but no 
general comparison seems possible short of numerous 
ex])eriments, 

3.2. Application to the Coverage Problem 

We observed in section 2.2 that the function 
F(X) we want to maxhnize could not be expressed 
in any neat formula. In fact, where jnore than two 
disks are involved, the construction of an algorithm 
to calculate the function would probably be too diffi- 
cult and time-consuming to be worthwhile. Wliat 
we need, first of all, is an approxunation for F(X). 
We assume that the circuLu* area A is centered at 
the origin and its radius R is a positive integer. 
Furtlu^rinore, the radius r of the disks is also a 
positive integer. ^^ The points of the plane (jfyO) 
where both p and q are integers are called (jrid- 
points and the area of A is appi'oximated by the 
number of grid-points Ui which lie inside the bound- 
ary of A, We shall also require that each pair 
{Xi,Yi) determining the center of a disk Ci be a 
grid-point. Our estimate for tlu* portion of A 
covered by the disks is the number of gi-id-points 7^2 
that lie inside A and at least one of th(* disks (see 
fig. 5). The ratio F{X) is approximated by tlie 
quotient 712/%. 

The problem is now to maximize F(X) where F 
is given by the approximation, over all vectors X 
such that the components are integers j in the 
range —R^j^ R. In this form the probleju satisfies 
conditions (1) and (2) of section 3.1. To satisfy the 
third condition we must specify the neighborhood 
N(X) for each vector X in the domain of F. We 
shall define N(X) to be all those vectors in the 
domain of /'which can be derived from J^by adding 
± 1 to exactly one component of X. This means 
that each Xhas 4?! neighbors except in the boundary 
situations (some component of X is ±R) where it 
has less. 

The problem as now formulated can be submitted 
to the ]nethods of maximization described above. 



11 Since the function FiX) is a ratio, the value of F is not changed if R,r, and 
the vector X are multiplied by a constant factor. 
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Some rather slight deviations from the general 
method were employed because of certain peculi- 
arities of our problem, but for the most part these 
methods were the techniques that were programed 
for use on a computer. 

3.3. Drawbacks of the Methods 

At this point we should acknowledge several 
rather sobering facts about our method. The black 
box maximization methods described (steepest, 
slowest, and first positive ascent) all find a point ^^ 
which is a relative or 'local'' maximum of the 
function; if one were certain beforehand that the 
function has only one such 'local'' maximum, there 
would be no further problem. Unfortunately, this 
is usually not the case because one is dealing primarily 
with functions whose behavior is generally unknown. 
The best way to increase the probability of hitting 
on the true ''global" maximum seems to be to 
repeat the search many times with random initial 
points and different schemes of ascent. Secondly, 
the discretizing of the problem which was effected 
in order to be able to approximate the function has 
introduced some error into the numerical results. 
Although theoretically the mesh can be refined ^^ 
to obtain any accuracy desired, the limitations set 
by time and the size of computer memory make it 
impossible to refine the mesh indefinitely. 

3.4. A Related Topic 

A topic related to the methods discussed in 
section 3 of this paper is that of maximizing an 
"unknown'^ function whose every evaluation requires 
physical experimentation and so, besides being 
costly, involves experimental errors. Since the 
classical paper of Box and Wilson [151 appeared, 
much work has been published in statistical journals 
on the design of efficient explorations schemes for 
such "response surfaces"; we mention here only 
a paper of Box and Hunter [16] and those by 
Brooks [17]. 

4. The Computer Program 

4.1. Specializations and Subroutines 

The technique we have described was programed ^^ 
and coded in FORTRAN and SAP for use on an 
IBM 704 electronic computer. The code has been 
debugged and a large body of data has been collected. 
There are several things which did not appear in the 
preceding description of the method but were neces- 
sary additions or at least were clearly indicated. 

The procedure for choosing an initial "vector" X 
to begin a search was left unspecified in the foregoing. 



12 " Point" is used here to mean an element in a vector space; i.e., a vector. 

13 The mesh is refined l^v multiplying the appropriate variables by a constant 
factor; this does not sound liVe refinement but it amounts to the same. 

1* See appendix 3; the FORTRAN program there reproduced is for the Cycling 
First Positive Gradient described in section 4.2, paragraph 5. 



As the computer code was written, the selection is 
made as follows : 

(1) Using "steepest ascent," the initial Xis chosen 
by randomly generating a certain number of vectors 
and selecting the one that achieves the highest func- 
tion value. 

(2) Using either "slowest ascent" or "first posi- 
tive ascent/' the initial X is selected by a single 
random generation. 

(3) Using any of the methods, the initial X may 
be read into the computer as an input variable. 

Another variation on the general method that was 
present in the program was refinement of the mesh. 
In section 3.2 we required that R and r be positive 
integers and that the centers of the disks be grid- 
points (i.e., Xi and Yi must be integers for all i). 
The value of R is further limited by the program to 
the powers of 2, and each time a relative maximum 
is achieved with some value of R the mesh is refined 
by doubling i?, r, and X. This is effectively the same 
as if we had actually refined the mesh of grid-points 
by constructing new lines halfway between those that 
already define our grid-points. After this refine- 
ment has been carried out, the search is continued 
until a relative maximum is found. The mesh is 
then refined once again and so on until the mesh is as 
fine as we desire. Since the fineness of the mesh is 
indicated by the value of i?, we specify as an input 
to the program the maximum value of R indicating 
the final mesh size. The reason for this successive 
refinement is that if the search is begun with a coarse 
mesh, the bulk of the searching process can be done 
in a lesser amount of time; this is because the amount 
of computation depends very strongly on the num- 
ber of grid-points, as might be expected. 

There are ten subroutines that are called for by the 
main program; a listing of these and a brief descrip- 
tion of each follows : 

MK2T — computes and stores in memory a table 
of the squares of all positive integers less than 1,000. 

XBAR — ^randomly generates n pairs of coordi- 
nates (X,, FO where X^+FK^' with X^ and 
Yi integers. 

SUMX — computes for the current X and mesh 
size R, the number of grid-points which lie inside 
the large circle A and at least one of the disks. This 
number is called NSUM. 

RAT — computes the total number of grid-points in 
A and divides NSUM by this number to get RATIO, 
our7^(X). 

VECTOR — for each particular value of r, selects 
those points near the boundary of a disk which 
should be scanned to determine whether some move- 
ment of a disk represents a gain or loss in NSUM.^'" 

NEAR — for each disk Cj a determination is made 
of those disks near enough to Cj so that they might 
overlap Cj if Cj were moved by one unit in some 
direction. For example, in figure 1 it should be 
quite clear that disk Ci need not be considered when 
we are interested in knowing what effect small 
changes in the position of Cr^ have on the value of 
NSUM. 



!■' See section 4,2. paragraph (4), for details. 
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SCAN 1— coinputes the change in NSUM due to 
one of the four possible changes in the position of a 
particular disk (), wlien the disk is entirely within 
the boundary of A and would still be so even after 
one of the four changes in its position. 

SCAN 2 — performs the same computation as 
SCAN 1 when the disk Cj is partially outside A or 
might be after a single move. 

XNEW — changes the 'Vector'' Xi to the new 
'Vector'' X2 indicated by the search procedure. 

KEFINE — nudtiplies the variables R, r, and X 
by two, thereby effecting the mesh refinement. 

For each of the different methods of ascent a 
separate program has been written. Each program 
uses the ten subroutines to do the bulk of the work. 

4.2. Shortcuts in the Program 

Some special techniques and shortcuts were used 
in the program and subroutines; they were devised 
partly to save time in computation and partly as 
a result of certain peculiarities of the specific problem 
to be solved. 

(1) The restrictions that were imposed on the 
problem in section 3.2 may seem somewhat artifi- 
cial. As a matter of fact, there was a definite rea- 
son for recasting the problem in such a way that 
almost all the variables involved are integers. The 
computer for which the program was written (IBM 
704) has separate sets of instructions to deal with 
integer variables and noninteger variables, and the 
time required to add two integers is two machine 
cycles ^^ where the (floathig) addition of two non- 
integers takes from 7 to 11 cycles. It was thought 
that this difference would affect greatly the total 
time required for a run. 

(2) Tlie computation of the table of integer squares 
up to 1000 which is done by the subroutine MK2T 
was inserted into the program to save time also. 
A single multiplication takes 20 machine cycles on 
the IBM 704 whereas looking up the square of an 
integer from a table stored in the memory takes 
only four cycles. 

(3) In the calculation of NSUM, tlie total num- 
ber of grid-points that are covered by the configura- 
tion X, it is possible to simplify the computation 
procedure by the following shortcut: 

The calculation of NSUM is performed by scan- 
ning all the grid-points on some vertical line x=k, 
an integer. The lowest grid-point on the line that 
is also inside A is the first to be scanned. For each 
successive grid-point proceeding in the direction of 
positive y values, a decision is made as to whether 
or not the point falls in one of the disks Cj. Sup- 
pose that a certain grid-point (X,Y) is determined 
to be inside a disk Cj with center (Xj,Yj), as in 
figure 6. Rather than contuming on to the next 
adjacent point, we calculate tiie quantity (Yj—Y) 
and draw the conclusion that all points between 
(X,r) and 

(X,Y+2{Yj-Y)) = (X,2Yj-Y) 



16 A machine cycle in tlie IBM 704 requires 12 ^sec. 



are also hiside Cj. The numbei- of points to be 
scanned has thus been substantially reduced. 

(4) The subroutine NEAR (see sec. 4.1) saves 
some time in computation by selecting from the 
set of disks, those that have no effect on small 
movements of a certain disk Cj. For example, if 
the program were about to determine what changes 
occur in NSUM as a result of moving disk Cq of 
figure 1 to any of the four positions possible, it 
would not be necessary to consider disks 6\, C2, C-iy 
or 6\ in making this calculation. WJien the total 
number of disks is larger this shortcut in the 
computation should be quite effect iv(\ 

(5) The general method calls for computhig the 
function value F{X) for all vectors X that are 
'^neighbors'' of the current vector Xi. From our 
definition of neighbor we can see that a neighbor 
of Xi corresponds to a configuration derived from 
that of Xi by moving some Cj one mesh-unit in one 
of four directions. Fortunately it is not necessary 
to compute F{X) at all neighboring points. Wliich- 
ever metliod of ascent is to be used, the important 
value to be computed is AF=F(X)—F(Xi), and for 
each X this can be computed without computing 
either F(X) or F(Xi). We ask the questions: How 
many mesh-points that were not covered by Xi are 
covered by X? How many points that were cov- 
ered by Xi are not covered by X? The answer to 
the first question tells us how many points liave 
been gained and the answer to the second, liow many 
have been lost. The difference between these two 
is the net gain in covered points due to changing 
Xi to X. In figure 7 the two separate sets of mesh- 
points indicate respectively the set of points that 
could either be 'lost'' or ^'gained" by moving the 
disk to the right. Furthermore, for any direction, 
two similar sets of mesh-points can be selected to 
scan in computing AF. 

(6) Using the terminology of the preceding section 
(5), only points that are within A are candidates for 
classification as gains or losses. Therefore, if the 
disk Cj being scanned is inside A by at least a mesh- 
unit (that is, all points of Cj are inside A and no 
closer than one mesh-unit from the boundary of A)^ 
then all the points on the periphery of Cj are inside 
A and qualify as possible gains or losses. In this 
case SCAN 1 is used to compute AF. However if 
Cj overlaps A, we must determine for each mesh- 
point being scanned whether it lies inside A. In 
this case SCAN 2 is used. 

4.3. Inputs and Outputs 

Very few input variables are required for the pro- 
gram. A list of the most important and a description 
of each follows :^^ 

N7 — the number of separate cases to be done. 

J7 — ^this variable indicates whether or not the 
initial 'Vector" is to be read in as an input. 

KO — Initial radius of large area A (this must be 
some power of 2). 



17 Integer variables in a FORTRAN code must be designated by a symbol 
beginning with one of the letters I, J, K, L, M, N. This is why the radius R 
of the disk A to be covered is denoted by KO. 
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NC — The number of covering disks. 

N — The number of random selections of a ^Vector'' 
to choose the initial vector. 

KBIG — This number represents the finest mesh 
to be used. 

JXPO— Satisfies the equation 2'^^o=KO. 

IKO — Initial radius of covering disks (satisfies 
IRO<KO). 

N4 — Number of times each case is to be repeated 
with different initial ^ Vectors.'^ 

As far as outputs are concerned, a fairly readable 
format has been devised. Initially, a general 
description of the case to be done and the important 
inputs are printed out. Next the initial 'Vector'' 
is printed out with its function value RATIO. 
Each move that occurs is printed out and the final 
configuration at each mesh size is recorded with its 
RATIO. The total number of moves that have 
been made is also printed out when the final con- 
figuration has been attained. 

5, Results of Computer Runs 
5. 1 . Description of Cases Studied 

A systematic series of runs was completed of 
cases involving between 2 and 10 disks where the 
ratio rjR took on interesting values^^ between and 1 . 
Table 1 gives the fraction of A that was covered 
by the best configuration found during the search for 
the global maximum. Diagrams of these con- 
figurations can be found in figures 8 to 44. These 
numbers are accurate to about ±0.002. A detailed 
explanation of the convergence properties of the 
approximation is given in section 6. 

5.2. Closeness of Relative Maxima 

In the cases where two, three, or four disks were 
used to cover, the final configurations were all global 
maxima. In the other cases there were as many as 



18 A case is not interesting if complete coverage is possible with a smaller number 
of disks, or if all disks can be placed inside A so as not to overlap. 




Figure 8. Case {2,9/16). 



five different local maxima found. Two configura- 
tions were considered to be different if the con- 
figurations were geometrically unlike. For example, 
the configurations of figures 45 and 46 are considered 
to be the same whereas both are different from that 
of figure 47. 

Diagrams of the local maxima that were found in 
several cases are depicted in figures 45 to 60. The 
value of RATIO is included with each configuration. 
It should be pointed out that the three best con- 

Table 0. Mean coverage ratio'^ with case (n,r) 





2 


3 


4 


5 


6 


7 


8 


9 


10 


Me 














0.373 
.461 
.541 
.610 


0.409 
.501 
.583 


0.442 


% 








0.321 

.385 
.445 
.500 
.551 


0.371 
.442 
.507 
.565 


0.418 
.494 
.562 


.538 


Me 






0.322 
.376 
.426 
.473 
.516 
.556 


.622 


|/2 


"a 242' 
.274 
.304 
.334 
.362 
.388 
.414 
.437 


0.298 
.341 
.381 
.420 
.456 
.490 
.521 




^16 






% 










11/16 












¥i 














iMe 














% 
















iMe 
















1 





































*The centers of the disks are randomly chosen from a uniform distribution on 
the circular disk of radius (1+r). 

The entry in the table is the expected value of the ratio of coverage where n 
disks of radius r are randomly placed as described. 

This table contains the mean coverage for the same cases for which table 1 gives 
the maximum coverage. 



Table 1. Maximum ratio of coverage with case (n,r) 
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1.000 
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1.000 


1.000 
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1.000 


X 


X 


X 


X 


% 


.686 


.883 


.982 


1.000 


X 


X 


X 


X 


X 


Hie 


.762 


.935 


.999 


X 


X 


X 


X 


X 


X 


% 


.829 


.972 


1.000 


X 


X 


X 


X 


X 


X 


n^(^ 


.889 


.994 


X 


X 


X 


X 


X 


X 


X 


Zi 


.939 
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1.000 
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X 
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— indicates all disks can be packed into large circle. 
X indicates total coverage is possible. 




Figure 9. Case (2,618). 
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Figure 10. Case (2,11/16). 




Figure 11. Case (2,314). 




Figure 12. Case (2,13/16). 
657118—62 3 




Figure 13. Case (2,7/8). 




Figure 14. Case (2,15/16). 




Figure 15. Case (3,1/2). 
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Figure 16. Case {3,9/16). 




Figure 17. Case (3,5/8). 




Figure 18. Case (3,11/16). 




Figure 19. Case (3,3/4). 




Figure 20. Case (3,13/16). 



figurations for the case ^^ (6,1/2), figures 48 and 49, 
differ in function values by at most 0.006. In case 
(10,3/8), figures 43 and 44, the two best local maxima 
differ by less than 0.002, and in case (9,5/16), figures 
57 and 58, with the mesh refined to 256, the two 
local maxima differ by less than 0.0002. In the 
latter case the ratios must be considered as in- 
distinguishable due to the error of the approximation 
at this mesli.^^ 



19 The following notation specifies the situation n=Q, r=l/2 (and R = l). This 
method of denoting cases will be used consistently in the rest of the paper. 
23 See section 6. 
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FiGUEE 21. Case {3,7/8). 




Figure 22. Case UJ/l^). 

In general it was found that the differences be- 
tween the local maxima in each case were on the 
order of one or two percent of the total coverage. 
The largest such difference that was observed was 
in case (8,%) where peaks of height 0.905 and 0.947 
were found, constituting a difference of about 4 
percent. 

5.3. Details for the Six-Disk Case 



A more detailed study was made of the local 
maxima achieved with six covering disks at radii of 
/i6, Vs, ^6, }2, and %6. The case (Q,%q) is one that we 
termed *'not interesting'' in section 5.1. because 




Figure 23. Case {4,1/ 




Figure 24. Case (4,9/16). 




Figure 25. Case {4,5/8). 
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Figure 28, Case {6,7116). 




Figure 31. Case (6,3/8). 
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Figure 32. Case (6,7/16). 




Figure 33. Case {6,1/2). 




Figure 34. Case (7,3/8). 




Figure 35. Case (7,7/16). 




Figure 36. Case (8,5/16), 




Figure 37. Case (8,^ 
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Figure 38. Case (8,7/16). 




Figure 39. Case (9,5/16) Ratio = 0.82U, 




Figure 40. Case (9,5/16) Ratio = 0.8241. 




Figure 41 . Case (9,3/8) . 




Figure 42. Case (10,5/16). 




Figure 43. Case (10,3/8) Ratio = 0.9924. 
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Figure 45. Case (10,3/8) Ratw = 0.99L 




Figure 46. Case (10,3/8) Raiio = 0.988. 




Figure 47. Case (10,3/8) Rafio = 0.990, 




Figure 48. Case (6,1/2), 




Figure 49. Case (6,1/2), 



195 




Figure 50. Case {6,112). 




Figure 51. Case {6,112). 



all of the disks can be placed inside A so as not to 
overlap each other; thus the ratio of coverage can be 
calculated accurately. The case (6,%6) is also ''not 
interesting'' because total coverage is possible. 
Figures 48 to 52 and 61 to 67 contain diagrams of the 
local maximum configurations for these cases. 

An interesting phenomenon occurs in these cases. 
The best placement of the disks in case (6J{6) is a 
''central'' or "flower-petal" configuration (see fig. 64), 
but in the case (6,K) the same "central" configuration 
is merely the third best. This seems to indicate that 
there is an intermediate value of r between Ke and K, 
where the "central" and "triangular" (see fig. 48) 
configurations both cover equal portions of the total 




Figure 52. Case {6,1/2). 




Figure 53. Case {10,3/8). 



area. The "ring" (see fig. 63 or 65) configuration of 
case {Q,yu) seems to be a fairly natiu*al analogy to 
either the "triangular" configuration or the "dia- 
mond" (see fig. 51) configuration of case (6,)0. 

A further comparison was made between the 
"triangular" and "diamond" configurations. Three 
cases were used— (6/%2=0.531); (6,6^(28=0.539); 
(6,^% 4= 0.547). In each case the "triangular" con- 
figuration of disks was the best, but its margin of 
victory decreased as the value of r increased. Table 
2 contains the results of the comparison.-^ 



21 The final mesli was 256 for all cases. 
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Table 2. Comparuou of two 6-disk con figurations 



r 


Triangular 


Diamond 


Difference 


0.531 
.539 
.547 


0. 99591 
. 99802 
. 99932 


0. 99496 
. 99768 
. 99928 


0. 00095 
. 00034 
. 00004 



Table 3. Comparison of two 10-disk configurations 



r 


Central 


Other 


Difference 


0.375 


0. 9924 


0. 9911 


0. 0013 


.383 


.9973 


.9957 


.0016 


.391 


.9996 


.9985 


.0011 


.398 


1.0000 


.9998 


.0002 


.406 


1. 0000 


1. 0000 


.0000 




Figure 54. Case (10,3/8), 




Figure 55. Case (10,3/8). 




Figure 56. Case (10,3/8), 




Figure 57. Case (9,5/16). 




Figure 58. Case (9,5/16). 



197 




Figure 59. Case (8,S 




Figure 60. Case (8,3/8). 




Figure 61. Case (6,5/16). 




Figure 62. Case (6,3/8). 



(6,'/e) 




Figure 63. Case (6,3/8). 



The chart shows that both configurations achieve 
total coverage at approximately equal values of r. 

Runs were also made of several cases involving 
10 disks at radii between %=0.375 and i%2=0.406. 
The two types of relative maxima were compared 
at five different values of the radius r. Diagrams of 
these relative maximum configurations for r=% are 
contained in figures 45 and 47. In all five cases 
the ^^centraP' configuration (fig. 47) was the best 
placement but usually by only a tenth of a percent. 
This data appears in table 3. 

Once again the two types of relative maxima 
achieve total coverage almost simultaneously. 

5.4. Comparison With Analytical Solution 

A run was made to compare our method with the 
analytical results of Neville for the case of five disks 
which was mentioned in section 2.3. According to 
Neville, the smallest value ot r for which coverage 
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Figure 64. Case (6,7/16). 




Figure 65. Case {6,7/16). 

with five disks is possible is 0.609375; tliis is called 
the ''critical radius.'' The results of the run were 
quite suiprising. The ratio of coverage was 0.999942, 
just less than total coverage. It had not been 
suspected that the method we used would achieve 
a configuration exactly like that of Neville (see 
fig. 3), since other runs seemed to show that the 
ratio of coverage for configurations that are associated 
with almost total coverage is rather insensitive to 
changes in the configuration. This made it seem 
likely that the configuration achieved by the com- 
puter run would not resemble the Neville configura- 
tion too closely. The configuration arrived at by 
our searching teclmique (see fig. 68) bears a remark- 
able resemblance to that of Neville. Each con- 
figuration contains three disks that intersect in a 
point very near the center of A and two others 
symmetrically placed with centers substantially 
displaced from the center of A. 

An analytical solution to the covering problem 
for the cases where n=2 and }i<Cr<^l is contained 




Figure 66. Case (6,9/16). 




Figure 67. Case (6,9/16). 
Note: Figure 68 is located on p. 182. 

in appendix 2. A comparison of the numbers gotten 
from computer runs with the true values showed 
excellent agreement, usually differing by less than 
0.002. 

5.5 Comparison of Search Methods 

A comparison was made of different methods of 
search to determine the Irequency of occurrence of 
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FiGUKE 69. Neighborhood of a mesh point. 




FiGUKE 70. Maximum "packing. 




Figure 71. Wall interrupted coverage by three disks. 




Figure 72. 1 disk coverage of wall region. 




Figure 73. 2 disk coverage of vmll region. 




Figure 74. 3 disk coverage of wall region. 




Figure 75. 4 disk coverage of wall region. 
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'^globar' maxiniji and to <>;el sojiu^ i(l(\i of the liino 
required by the tlifrereiil methods. Five ''iiu^thods" 
were prograiued and run. Tliey are all derived from 
the three deserihed in section 3.1, witli minor 
variations. 

(1) Plateau Steepest Ascent : This is basically the 
^'steepest-ascent" method with two exceptions. 
First — ^if at some stage of the search there are no 
positive moves ^^ indicated for the current configu- 
ration X, then the neighborhood of X is effectively 
enlarged hy allowing each disk to move to any one 
ol the eight grid-pohits surrounding its center (points 
A through H in fig. 69). This is done on\y wlien 
no moves are possible under the original definition 
oi '^neighbor.'' Second — if both of these searches 
fail to obtain a positive move then a search is made 
to determine if all four moves (that is, a move from 
to one of B, D, F, H in fig. 69) of some disk are 
zero moves (i.e., leave RATIO unchanged). If 
such a disk is found, successive random placements 
of its center are tried until the value of RATIO is 
increased. If 25 trials fail to accomplish anything, 
the searcli is abandoned, and the refinement of mesh 
is made, etc. 

This second feature is what gives the pi'ocedure 
the name 'Tlateau." In a function of a shigle 
variable yC.r) one says that a point .x,, is on a plateau 
of the function if all points in some neighl)orhood of 
Xq have equal function values. In two dimensions 
a plateau, analogously defined, can be visualized as 
being a flat portion of the surlace represented by the 
function (j{x,y). An intuitive idea of the meaning 
of ' 'plateau" is possible in these two cases precisely 
because a ^'picture" of the function can be repre- 
sented in three dimensions or less. 

Although we cannot visualize ])lateaus in Jiigher 
dimensional spaces, we nevertlieless defme plateau 
analogously as a neighborhood of t he space on whicJi 
the function is defined throughout which the 
function is constant. In our case we have considered 
F as a function of two variables, tlie (coordinates of 
the center of a single disc, and so we have the type 
of plateau that makes sense visually. 

It should be mentioned here that this particular 
feature of the 'Tlateau Steepest Ascent" method 
did not have nearly the same significance as did 
the first feature (the enlargement of neigliborhood). 
That is to say, the results in using the steepest 
ascent method were more drastically altered by 
neighborhood enlargement than by the ' 'plateau" 
feature. 

(2) Least Positive Ascent: This is the ''slowest 
ascent" described in section 3.1. 

(3) First Positive Ascent: As described in section 
3.1. 

(4) Cycle First Positive Ascent: This is the same 
as (3) except that the first disk whose moves are to 
be tried is the disk immediately following the disk 
that was last moved. The disks are ordered in a 
cycle Ci, G, . . ., Cn-u Cn, Ci, etc., for this purpose. 
This variation of (3) was used because the return to 



22 A positive move from a configuration X is effected by changing X to one of 
its neighbors Xi such that /'XXi)>i^(X). It is accomplished by a change in the 
position of a single disk in one of the coordinate directions hy one mesh unit. 



Ci at (MicJi stage seemed to hitroduce some bias that 
was undesirable. 

(5) (Vch^, First Positive Ascent with Eight De- 
grees of Freedom: this resembles (4) except that 
tlu'ougliout the search the neighborhood of X is 
the enlarged neighborliood employed by the Plateau 
Steep Ascent. 

Two cases, (6,K) and (10, JO, were selected to be 
used in the comparison. Each case was known to 
have nonglobal local maxima that occurred quite 
frequently, and the configurations were quite distinct 
from a visual standpoint. Table 4 describes tJie 
frequency of occurrence ^^ of the global maximum 
in both cases, according to which method was used. 
The last column is probably the most significant 
because the cost depends on the time consumed and 
not the number of trials. 



Table 4 


Comparison 


of the 


5 ascent methods for 


2 cases 


Method 


Case 


Trials 


Global 
max- 
ima 


Time 


Global 
maxima 
per trial 


Time 
per 
trial 


Global 
maxima 
per min 


1- 


(6,^) 
(6,H) 

(6,H) 

(10,^) 
(10,H) 

(io,H) 

(10,3/8) 

(10,^) 


25 
25 
25 
25 
25 
25 
21 
24 
25 
20 


13 
10 
11 

7 
8 
5 


5 
5 


Min. 
44 
47 
34 
19 
52 
82 
89 
71 
40 
90 


0.52 
.40 
.44 
.28 
.32 
.20 
.00 
.00 
.20 
.25 


Min. 
1.76 
1.88 
1.36 
.76 
2.08 
3.28 
4.24 
2. 96 
1.60 
4.50 


0.295 




.213 


3 

4 

5 


.324 
.368 
.154 


1 


.061 


2 


.000 


3 


.000 


4 


.125 


5- . 


.056 







From this point of view, the Cycle First Positive 
Ascent (4) is clearly the best. There is another 
point of view, however, that may be still more 
significant in a comparison of these five methods. 
The point is this— the method to be preferred is the 
method that tends to reach more different peaks. 
Such a method would reach the global maximum 
less often but one would feel more confKhMit of the 
liighest peak among 5 or 10 than he would if only 2 
or 3 distinct peaks had been found. Table 5 con- 
tains a record of the number of distinct peaks that 
were found by each method. It seems to indicate 
that the variations of the ^'first positive ascent'^ 
method are better than the others at locating different 
peaks. 

Table 5. Number of different peaks found by the 5 methods 
Case (6,1/2) 



Method 

Peak's ... 


1 
3 

1 
3 


2 
4 

2 
3 


3 
4 

3 

4 


4 
3 

4 
5 


5 
4 


Method 

Peaks 


5 
3 







5.6. Critical Values of Radius 

Considering problems (c) and (d) of section 1.2 
once again, we define the ^low critical radius'^ of n 
disks ri{n) as the answer to problem (d) and the 



23 Two final configurations are deemed equivalent, in this comparison, if they 
are similar geometrically and achieve nearly identical values of RATIO. 



201 



'^high critical radius/' r2(n) as the answer to problem 
(c). Neville [1] refers to the latter simply as the 
^'critical radius'' but we want to study both. For 
any number of disks n the cases that are ' 'interesting" 
are the cases {n,r) where ri(n)<r<r2{n). 

A determination was made of ri(n) and r2{n) for 
all n in the range l<n<10 (table 6). Some of the 
values were found by using the search methods and 
others could be determined analytically. For ex- 
ample we can calculate ri(n) for 2<n<5 by the 
following argument: 

We want to place the n disks in a ring around the 
center of A in such a way that the disks are packed 
in as tightly as possible without overlapping. Under 
these conditions each disk requires a sector cut off 
by an angle d=27r/n. Referring to the diagram in 
figure 70 we can immediately write 

r/(l—r) = sin 6/2 and ^=27r/n, 



which reduces finally to 



sin (tt/ti) 
1+sin (ir/n) 



Table 6. Low and high critical radii for up to 10 disks 



n 


1 


2 


3 


4 


5 


6 


7 


8 


9 


10 






ri{n) 

r2{n) _ 

Ar 


1.0 
1.0 



0.50 
1.0 
.500 


0.464 
.866 
.402 


0.414 
.707 
. 293 


0.370 
.609 
.239 


0.333 
.555 
.222 


0.333 
.500 
.167 


0.302 
.437 
.135 


0.276 
.422 
.146 


0.266 
.398 
.132 











and this formula is valid so long as the ring configura- 
tion is clearly the optimum packing. 

For small values of n, T2{n) can be calculated by 
noting that to achieve total coverage, all of the cir- 
cumference of A must be covered. For a given n 
the minimum r required to cover the circumference 
is given by: 

r=sin tt/ti. 

If total coverage is in fact achieved at this value 
of r then r2{n) has been found. This argument was 
found to be valid for 7i=2,3,4. 

Table 6 shows that the length Ar of the interval 
(ri{n) ,r2(n)) shows a tendency to decrease with n, 
although not monotonicallj^ Thus for larger n 
there is a smaller range of interesting cases. In 
fact when n=10 the interval is (0.266,0.398) whose 
length is only 0.132, whereas when n=2 the interval 
is (0.50,1.0) of length 0.500. 

5.7. Efficiency of Covering Configuratioii 

It was thought that some measure of the efficiency 
of a covering might be useful in some applications. 
The efficiency should give some indication of ratio of 
coverage versus total covering area available. We 
therefore define the efficiency E{n,r) as the ratio 
between the maxunum area (not percent!) coverable 
by n disks of radius r, and the total composite area 



of the n disks. If we denote the best percent cover- 
age by C(n;r) and remember that we have been 
assuming R=l we obtain the formula for efficiency: 



E{n,r) = 



TrC(n,r) _jO{n,r) 



mrr^ 



nr^ 



The values of C{n,r) are contained in table 1, and 
E{n,r) in table 7. For most pairs {n,r) the efficiency 
is greater than K but two cases were found where it 
was less than K- Specifically, £'(3,7/8) = 0.435 and 
£(4,3/4) =0.444. 





Table 7. 


Efficiency of coverage, E(n,r) 






\ 




















r \ 




















2 


3 


4 


5 


6 


7 


8 


9 


10 


Me 


1 


1 


1 


1 


1 


1 


0.985 


0.934 


0.900 


% 


1 


1 


1 


0.989 


0.936 


0.907 


.842 


.772 


.706 


Me 


1 


1 


0.974 


.883 


.796 


.731 


.653 


.580 


.522 


Vi 


1 


0. 961 


.863 


.749 


.653 


.571 


.500 


X 


X 


91 e 


0.948 


.858 


.740 


.625 


.527 


X 


X 


X 


X 


% 


.878 


.753 


.629 


.512 


X 


X 


X 


X 


X 


iMe 


.807 


.659 


.528 


X 


X 


X 


X 


X 


X 


% 


.737 


.576 


.444 


X 


X 


X 


X 


X 


X 


i^e 


.674 


.502 


X 


X 


X 


X 


X 


X 


X 


% 


.614 


.435 


X 


X 


X 


X 


X 


X 


X 


1^6 


.556 


X 


X 


X 


X 


X 


X 


X 


X 


1 


.500 


X 


X 


X 


X 


X 


X 


X 


X 



X indicates total coverage is possible, so increases in P=nr^ are just wasted. 



Table 8. 



Comparison of covering efficiency with nr^ held 
constant 



Ratio.. 

wr2 

E{n, r) 



(4,3%4) (5, 

0.973409 0. 

1.485 1.495 

0.655 0.656 



%4) 



(6, ¥z) 
0. 979418 
1.50O 
0.653 



(7, 15/32) 

0. 995094 

1.538 

0.647 



(8,^6) 
0. 999432 
1.531 
0.653 



(9, i%2) 
0. 999130 
L485 
0.673 



(10, 2%4) 
0. 999626 
L526 
0.655 



Several runs were made to obtain data pertaining 
to problem (e) of section 2.2. The value of P=nr'^ 
was set at 1.5 (or as nearly as possible)^* and the 
corresponding efficiencies were calculated. At first 
it was suspected that efficiency would always be 
better when n was larger, but the results in table 8 
contradict the conjecture. The efficiency of (5,35/64) 
is 0.656 and that of (6,1/2) is only 0.653. This 
discrepancy could hardly be due to the errors in- 
volved in the approximation because the approxima- 
tion tends almost always to be an under estimate of 
the true value (see sec. 6). 

The Kerschner result cited in section 1.3 can be 
rewritten slightly to read: 



A 



lim ,^, . 9- 



9/27rV3=0.827. 



The expression on the left is the limit of the 
efficiency in the case of total coverage as the number 
of circles increases and r decreases. 

Denoting X= 2^/3/9 we get the following corollary 
to the Kerschner theorem : 



24 The possible values of r are restricted to numbers of the form /c/64 where A; is a 
positive integer <64. 
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Given a set S of area P/\ and e>0, tliere exists a 
number ro>0 such that |rW(r)— P|<^e when r^C^o. 

This follows immediately by writiiio- Kerschner^s 
theorem as: 

limrW(r) = X(P/X)=P. 

r->0 

If we define the potential P of a set of A^ disks of 
radius r as P=Nr^, the above corollary says that 
any set of area P/X can be covered by a set of disks 
with potential arbitrarily close to P. 

In particular, if the set is a disk of radius R, then 
we have 7rK^=P/\ or R= (P/TrX)^/' where l/7rX=-0.827. 
If we put P=1.5 we get P= 1.114 and this means 
that a disk of radius 1.114 can be ^^ilmost^^ covered 
by disks with a potential of 1.5. 

6. Convergence of the Approximation 

Some estimates of the accuracy of the approxima- 
tion used in the search methods have been 
determined. The approximation was introduced in 
the calculation of tlie ratio of coverage (RATIO). 
The convergence of the approximation as the mesh is 
refined was studied by an auxiliary computer pro- 
gram called RATIO CONVERGENCE which cal- 
culates the value of RATIO for a particular con- 
figuration at mesh sizes S, 16, 32, 64, 128. Several 
of the configurations were simple enough so that the 
ratio of coverage could also be cal emulated accurately 
by analytical methods. Table 9 gives the results of 
several runs perfonned to estimate the convergence. 

The configuration of the first case in table 9 con- 
sists of a single (hsk of radius r=l/2 centered at the 
center of A. This means that the ratio of coverage 
should be exactly 0.250. The table indicates that at 
a mesh of 64 the approximation differs by only 
0.0006 from the correct ratio of coverage. ^^ At a 
mesh of 256 the difference is only 0.00016. 

In the first three cases tiie change in RATIO 
between mesh 64 and mesh 128 was less than or equal 
to 0.0004. 

In the fourth case the configuration involved two 
disks and the correct ratio of coverage was computed 
using the algorithm referred to in section 2.2. The 
true ratio of coverage was found to be 0.40225, which 
means the approximation differs by less than 0.001 
from the correct value at mesh 64. 

Two other configurations were tried in wiiicli all 
disks lie entirely inside A and no overlap occurs. 
The cases were (6,5/16) and (6,2/8). In the former 
case the approximate ratio at mesh 256 was 0.584940 
and the true ratio was 0.5S5937, a difference of 
0.001. In the latter case the approximate ratio 
was 0.374503 and the exact value 0.875000, a 
difference of 0.0005. 

These tests of the accuracy of the approximation 
seem to indicate that at a mesh oF 64 one can expect 



Table 9. Comparison of ratio convergence at different meshes 



\ 

Case 


. Mesh 

\ 


8 


10 


32 


64 


128 


1 


0.233161 
. 865285 
. 860104 
. 378238 


0. 243380 
. 894073 

. 882724 
. 392182 


0. 247426 
. 897972 
. 885179 
. 398440 


0. 249436 
. 899681 
. 886295 
. 401354 


0. 249840 


2 _ _ 


. 899998 


3 


. 885882 


4 


. 402049 







25 A final mesh of 64 was the one used in most of the computer maximization 
experiments. 



accuracy ^^ on the order of ±0.002. Furthermore, 
the approximation was usually an underestimate 
of the true value. At a mesh of 256 the accuracy 
is usually about twice as good. Most runs used a 
final mesh of 64 since the additional time required 
to increase to a mesli of 256 did not seem to be 
warranted by the additional accuracy obtained. 

An attempt was made to obtain rigorous limits 
for the error involved in approximating the covered 
area by the number of grid-points. 

Let i(r), denote the number of lattice points 
inside the disk defined by x^ + if<t'', and let 

D(:r) = \Trr^-L{r)\. 

It is easily shown [18] that T){r) converges to 
zero as r-^oo^ but the question of just how fast 
is another matter entirely. 

Hilbert and Cohen-Vossen show that 
or in other words 

Landau [19, pp. 183-278] shows that 

/)(/') ^0(^2/3), 

in fact [19, p. 271], 

Z>(r)-O(r^^/^^+0 lor any €>0, 
but that 

Z>(r")^0(ri^2). 

These have been somewhat improved [20] to read, 
T){r)={){r''''^^^) for any e>0. 

The best result to date seems to be 

7)(r)=0(ri3/20) 

which was shown by Loo-Keng Hua in 1940 [21]. 
The ^'conjectured^' result is 

Z)(r)-:0(r^/2+.) for e>0 



26 This refers to the error in estimating F{X) by the approximation mlm defined 
in section 32. The error involved in estimating a relative maximum of F{X) by 
one of the final configurations produced by the computer program is somewhat 
larger. 
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and it has already been pointed out above that this 
is the best result possible. 

Obtaining numerical upper bounds for D{r) from 
these results would involve a great deal more effort 
than seems worthwhile here. The problem is clearly 
tied in with fairly abstruse number theoretical in- 
vestigations. Furthermore, the error bounds that 
may be derived from the Hilbert and Cohen-Vossen 
inequality are so bad that there is some real doubt 
as to whether the subsequent tightening of this 
inequality is substantial enough to help us out. 

There is also the problem of estimating the error 
for n possibly (and indeed probably) overlapping 
disks. The above results pertain to the estimation 
of the area of a single disk without any bites taken 
out of it, whereas we have a considerably more in- 
volved situation, especially since we really want as 
tight an inequality as possible. 

A recent paper by H. L. Mitchell III [22] gives a 
large amount of numerical results concerning L(r), 
the number of lattice-points (grid-points) in a circle 
of radius r. The paper includes calculations of D{r) 
and also D{r)/r^^^. The latter values were calculated 
to get some evidence for the conjecture that 
2)(r)=0(r^/-+0 for every e>0 mentioned above. 



7. A Conjecture Refuted by the Study 

During the research on the covering problem the 
following conjecture, essentially a ^'law of diminishing 
returns,^ ^ was formulated: 

If C(n,r) denotes the maximum ratio of coverage 
attainable with n disks of radius r, then 

C(n+l,r) - C(n,,r) < C(:n,r) - C(n-l,r), 

This means that the successive gains in coverage 
by the addition of disks one at a time are monotoni- 
cally nonin creasing. Unfortunately two of the cases 
that were studied produced results that contradict 
the conjecture. Table 10 gives the values of the 
ratio of coverage for the cases involved and the differ- 
ences. According to the conjecture the second dif- 
ference should be a negative number, but table 10 
shows the two cases found in which it is positive. 
The amount by which it is positive is large enough 
so that the error of the approximation could not be 
responsible for the sign of the second difference. 



Table 10. Results of several cases showing positive second 
difference in C(n,r) 



Case 



(5,3/8).. 
(6,3/8).- 
(7,3/8)-- 

(8,5/16). 
(9,5/16). 
(10,5/16) 



Ratio (mesh 
of 256) 



0. 700365 
. 794168 
. 895913 

. 773017 
. 824466 



First differ- 
rence 



+0. 093803 
+. 101745 



+. 051449 
+. 056352 



Second diffe- 
rence 



+0. 007942 



+. 004903 



8. The Second Covering Problem 

8.1. Description of Problem 

The second covering problem that was studied is 
somewhat more complex. We are given a circular 
area A of radius R and a certain number of straight 
lines which intersect the area A and divide it into m 
regions Ri(l<i<m). In addition, we are given a 
certain number n of circular disks of radius r<^R. 
The problem is to find that placement of the n disks 
which ''covers'^ the largest amount of the area of A 
subject to the following restriction: 

A point is considered to be covered if and only if it lies 
inside some disk whose center lies within the same region 
as the point in question. 

For example, the area of A that is covered by the 
three disks of figure 71 is shaded. Note that the part 
of disk Ci that is in i?i is not covered because Ci is 
not centered in Ri, but rather in i?2. 

8.2. Analysis of Problem 

A configuration of n disks that maximizes the cov- 
erage will necessarily have a given number ni of disks 
centered in each region Ri. Furthermore, the place- 
ment of the 7^^ disks in Ri constitutes the best cover- 
age oi Rihy ni disks independently of what occurs in 
other regions. 

This observation enables us to separate the prob- 
lem into two parts and solve it as follows: 

(1) For each pair (i,j) subject to l<i<m and 
1 ^j^ ^? calculate the maximum area of region Ri 
that can be covered by j disks. Denote this area 

by Mi,j)- 

(2) Let $ be the family of m-dimensional vectors, 

V={ni,n2, . . ., Um) such that the nt are all non- 



m 



negative integers and y^ n.-.^ 

i=l 
define the sum, 



-n. 



Now for each Ve^ 



SiV)=j:A(i,n^); 

i = l 

the maximum coverage we seek is then given by 

;S,ax=Max {S{V): Ve^} . 

Portion (2) of the solution is a purely combina- 
torial problem whose solution depends only on the 
values of the entries in the matrix A{i,j). We shall 
return to this problem later. 

8.3. Computing A[ij) 

According to the above formulation, the first step 
in a solution to the problem is the calculation of the 
matrix A{i,j) for l<i<m, 1 <j<n. For a particular 
{i,j) this means finding the maximum area of a spe- 
cific region Ri that can be covered with j disks. If 
the region Ri were circular, then the calculation 
would be that of the problem described earlier (sec. 
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1). This strong" similarity sug<;'est(Ml \hc ])ossi})ility 
of using a similar method to (•alciilute A(iJ). As it 
turned out, oidy very minor cluingos wore necessary 
to transform the methods of solving the earher prol)- 
lem into methods which will compute A{i,j). The 
changes essentially amount to restricting the centers 
of disks to grid-points inside Ki, and not counting 
points outside Ri when computing the coverage of a 
particular configuration. So far as the computer 
programs were concerned, these changes were effected 
with a minimum of difficulty, considering tlie usual 
complications which arise in modifying computer 
codes. The basic reason for this was that the new 
problem differed from the former only in the region 
to be covered, and therefore many of the complexities 
of the program were unchanged. 

The program accepts as inputs certain parameters 
specifying the lines that, along with the boundary of 
the large circle, form the boundaries of each region 
Ri in question. A line is specified as a 'lower slope" 
or an 'dipper slope'^ at input time, according as the 
region Ri lies above or below the line. As is well 
known, any line in the plane (excepting vertical lines) 
can be written in the form 

y=px+q. 

We assumed further that our lines have rational 
slopes and ^/-intercepts. That is, p=Pilp2 and 
^=5i/^2 where Pi,P2,(lug[2 are integers with :P2>0,g2>0. 
This will be true whenever the line passes through 
at least two grid-points.^^ Our equation then can be 
written as, 

y=^Pixlp2 + qi/q2 
or 

(P2(i2) y=(piq.2)x+(qip2), 

which is of the form 

ay=bx+c, 

where a,b,c are integers and a=p2g[2^0. These three 
integers a,b,c are the input parameters that specify 
a line L. If L were further specified as a lower 
slope for the region Ri, then any point (x,y) which 
lies in Ri must satisfy 

y>px+q^={b/a)x+(c/a), 

that is, ay>bx+c. It should be noted that this in- 
equality can be tested by the computer using integer 
arithmetic. This is precisely why we required that 
p and q be rational. For any (x,y), only a finite 
number of tests need be made to decide if {x,y) lies 
in Ri. 

An example of the results of using this program to 
compute A(i,j) for a particular region Ri(;r=%) is 
given in figures 72 to 75. The only interesting cases 
are j= 1,2, 3, 4 since total coverage would certainly be 
possible with five disks of the same radius. Since 
the computer program is essentially the same as 
that used in the solution of the former covering prob- 

27 Suppose the line y=px-{-q goes through points Pi = (r,.s) and P2 = (t,u); then 
it can be shown easily that p={u—s)lit—T) and fi = s—pr. Therefore if r,fi,t,u 
are integers, then p and g are rational. 



lem, there seemed no need for an exhaustive examina- 
tion of cases. The particular difficulties in the sec- 
ond problem may in fact lie in the combinatorial 
question described in section 8.2. 

8.4. Combinatorial Aspect 

Returning to the second part of the problem 
described in section 8.2 we find that the combina- 
torial problem can be reformulated ^^ as an integer 
linear programing problem. We sliall use a^ instead 
of A(i,j) below, where "i^' indexes subregions includ- 
ing a fictitious zeroth subregion to absorb any unused 
disks. Let the matrix (Xij),0<i<M, 0<j<N he 
defined as follows: 

Xij=l if exactly j disks are allotted to the ith 
subregion and x^-^^O otherwise. We can then express 
tlie problem as follows: 

Maximize 2^•, ^ aijXij subject to constraints 



(a) Xij>0 

(b) f:xij=l 

j=0 



0<i<M, 0<j<N 
()<i<M 



(d) Xij an integer 0<'i<iV/, 0<j<A^. 

As such, the problem can be handled at least in 
principle by the methods developed by Gomory [23]. 
Perhaps an especially effective algorithm can be 
constructed for the special problem involved here. 
This remark is added because the general method of 
Gomory has been found to converge unacceptably 
slowly in some cases. 

9. Appendix 1. Algorithm and Analysis for 
Two-Disk Coverage Formula 

9.1. Description of Parameters 

Suppose given a circle A of radius f=l, and two 
other circles B and C, of radius r<^l, whose centers 
are at respective distances di and ck from the center of 
circle A. Let 6 be the angle between 

(1) the segment joining the center of B to that 
of A, and 

(2) the segment joining the center of C to that 
of A (see fig. 4) . 

We construct an algorithm to compute the area ^^ of 
the ^ ^coverage set'' M which is common to circle A 
and at least one of the two circles B and (7, 
M=An{B[JC). 

We may assume that the interiors of B and Cboth 
meet the interior of A (i.e., both Af]B and A(]C are 
nonempty). Further restrictions on parameters 
{r,d,di,d>) are as follows: 

(1) 0<r<l 

(2) c?i— r<l; 6^2— r<l since Ar\B and Af]C are 
nonempty. 

2s This integer prosramine; formulation was suggested by A. J. Goldman (NBS 
Operations Research Section). 
29 \X\ shall mean- the area of X; thus area of M= \M\, etc. 
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Figure 76. Area common to two 





Figure 77. Area common to three disks. 




Figure 78 Calculation of angles 5, 6, 62. 




Figure 79. Computation of u and v. 
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Figure 80. Area hounded hy four arcs. 

(3) 0<^<7r since replacing d by (2^ — 6) or ( — 6) 
does not change area. 

(4) di<d2] we may assume for convenience that 
B's center is no fm'ther from ^'s center tlian Cs 
center is. 

We shall create an additional parameter c (de- 
pending on 6,di,d2)y 

c={d:i+di-2did2 cos ey^'; 

then c is the length of the segment joining the center 
of B to that of C as can be shown easily from the law 
of cosines. 

Tliere follows an analysis, case by case, of the 
various configurations that require different treat- 
ment when calculating the coverage area in terms of 
the given parameters r, 6, di, d2, and the defined 
parameter c. 

9.2. Intersection of Two Disks 

We begin by calculating the area common to two 
disks of radii r and R respectively with distance of 
centers di<Cr-\-R. In case Q;<7r/2 (see fig. 76a) 
the area may be calculated as follows: ^^ 

Area=Sector {QiOiQ2)-AQiOiQ2 + Sector (Q2O2Q1) 
-AQi02Q2, 

Sector (QiOiQ2) = (2l3/27r)7rR'=R'p, 

Sector (Q202Qi) = {2a/2T)7rr^=r^a, 

AQiOiQ2+AQi02Q2 = 2(Arh02Q2)=diR sin p. 

We finally get the formula 

AYesi=R^P+r^a-dJi sin 13. 







30 The symbol AQ1O1Q2 means the triangle with these vertices. Sector (Q1O1Q2) 
shall denote the area of the sector taken in a clockwise sense from Qi to Q2 about 
Ou 





e f 

Figure 81. Geometry of two disk cases. 

In case a>7r/2 (see figure 76b) the area is: 

Area=:Sector (^202^1) + (Sector {Q1O1Q2) 

-2(AO/>2ft)). 

This formula reduces to the same formula as the 
case a<^7r/2 by a similar argument. 

Furthermore in case a = Tr/2 our formula gives 
Area=i?^iS+7r/'72 — f/ii? sin 13 which is correct also. 
We have thus shown that for 0<Ca<7r, 

AYe2i=R^^+7'^a-diR sin /3. (1) 

Use of the law of cosines on AO1O2Q2 yields 

r^=d\+R^-2d,R cos fi 

R^=dl+r^-2rd, cos a, 
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whence we get a and j5 as 

a=arc cos {{r^+d\-B^)l2rd,) 

P=^YC COS ((R'+dl-r')/2Rd{). 

In particular we notice that if both disks are of 
equal radius r, and di=c, then we get 



Area=2r^ arc cos (c/2r)—cr sin (arc cos (c/2r)). 
If R=l and r<^lj the common area becomes 
Aresi=^ + r^a—di sin 13, 
where a and are given as 

a^arc cos ((r^+c?? — l)/2rc/i), 
/5=arc cos (il+dl-r')/2di). 
9.3. Intersection of Three Disks 



(2) 

(3) 

(4) 
(5) 



We now calculate the area common to three disks 
where the common area is bounded by three circular 
arcs (see fig. 77a). In this figure Ai,Bi, and Ci 
represent the respective centers of A,B, and C and 
P, ^,5" represent the points of intersection of the three 
circular arcs. The common area will be calculated 
by adding together the area of APQS and the three 
areas each bounded by one of the arcs and its 
associated chord. We call these slivers (see fig. 77b). 

First we must calculate /3, the angle subtended at 
Ai by the arc PS of disk A. Using the law of cosines 
on AAiBiS and AAiCiP, remembering that BiS=r 
= CiP and putting ^^ iZ = 1, we get 

a+/3=arccos {{dl-r^+l)l2di) 

^+^=arccos {{di-r'+l)/2d2). 

Adding the two equations and subtracting the 
equation a-\-p-\-y = S, we get 

/5=arc cos ((df— r2H-l)/2rfi) + arc cos ((dl—r^ 

+ l)/2d2)-d. 

Kef erring to figure 77b, we calculate the area of 
the sliver associated with angle (3 as the difference 
between the sector (PAiS) and AAiPS. Denoting 
the area of the sliver as S^ we have, 

S0 = Sectoi {PA,S)-AA,PS 
= {Pl2ir)ir-xy 
= /5/2 — sin (/5/2) cos (^/2) 
= (/3-sin/3)/2. .,... (6) 

We must now calculate the angles b and e of 
figure 77a. We shall present the argument leading 

31 This assumption involves no loss of generality since the ratio of coverage oiA 
depends only on (r//?). 



to the calculation of b; the calculation of e is similar. 

Referring to figure 78a we wish to calculate the 
angle ^QC^P=b^ + b2 = d. Notice that P might lie 
on the left side of BiCi, say at P', in which case the 
angle ^2 if measured by (^4—53) would be negative. 
The formulas we shall derive will not be affected hy 
this difference, as can be readily checked. 

We note immediately that 

b = b,-\-b,-b, (7; 

so that we need only to calculate bi, bs, b^. 

Using the law of cosines or sines on the proper 
triangles and noticing that AQCiBi is isosceles we 
obtain 

cos bi = {c/2)/r=c/2r, 

sin bs/di=sin d/c, 

l^r'^j^(ll—2rd2 cos 54. 

Solving for bi, bs, 54 we get 

5i=arc cos (c/2r) (8) 

63= arc sin (di sin d/c) (9) 

54 = arc cos ( (r2 + (/2 _ 1 ) /2rd2) ( 1 0) 

Substituting into eq (7) we get 

5=arc cos(c/2r) + arc cos((r^+(i|— l)/2rc?2) 

— arc sin {di sin d/c). (11) 

We get a similar expression for e, 

e==arc cos((r^+(ii — l)/2r(ii) + a,rc cos(c/2r) 

— arc sin ((^2 sin d/c). (12) 

An argument similar to that leading to eq (6) can 
be used to show that 

S,=r'{b-sm b)/2, 
Se=r\e-sm e)/2. 

Referring to figures 77a,b we can readily see that 
the sides Si,S2,Ss of APQS are given by 

81=2 sin (i6/2) 
S2 = 2r sin (5/2) 
s^ = 2r sin (e/2). 

The area of APQS can then be gotten from the 
semip crime ter formula 



where 



(13) 



Area A='yjs(s—Si) (S—S2) (s—s^) 
s={si+S2+Ss)/2, 
In our case s=sin(/3/2)+r sin (5/2)+r sin (e/2), 
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and the coiniuon urea of the tliree disks is 

+ (l/2)(^-sin p+r^id-^in 5 + e-siii e)). (14) 

9.4. Additional Parameters 

In the algoritlini which is to follow shortly it is 
necessary to calculate the distance frojn Ai, the 
center of ^4, to the two points, Q and Qo, where the 
disks B and C intersect. We label these points so 
that Q is the closer to Ai and refer to their distances 
from Ai as u and v, 



u=AiQ; v=AiQo; u<v. 

If we consider the disks B and C as already 
occupying; fixed positions in the plane and consider 
all possible positions of Ai we notice almost innne- 
diately (referring to fig. 79) that Ai must lie below 
BiCi and to the left of QQq. This is because di<d2 
and u<v. We distinguisli two cases. 

If Ai lies in the region Ri, we have tlie following 
equations: 

u^=r'+di-2rd2 cos <p2, (0<(^2<7r/2) (15) 

<^i + «p2 = arc cos (c/2/0, (0<(^o<7r/2) (16) 

sin (^1/^/9 = sin d/c (17) 

and these can be solved for tlie value of u~ as: 

u^ = r^-\-(l^—2rdo cos [arc cos (c/2r) 

— arc sin (r/i sin d/c)]. (18) 

If A] lies in region Ro, say at A[ in the figure, then 
we get the equations: 

u^=r'+dl-2rd2 cos <P2, (0 <<^2<7r/2) 

(Pi — iP2 = ^YC COS (c/2r), (0<^i<7r/2) 

sin (pi/di=sm d/c 

and these can be solved for the value of u- as: 

^2^^.2_|_ j2_2^^2 COS [arc sin (di sin d/c) 

— arc cos (c/2r)]. (19) 

We notice that the expressions in brackets in eqs 
(18) and (19) differ only in sign and since cos x is an 
even function (cos (— x)=cos x) the formulas for u^ 
are identical. 

In the first case (^1 in Ri) we liave the following 
expressions 



^2^^2_|_^2_2^^2 COS ((Po+(Pi) 

where ipo is given b}^ 

<^o=arc cos {c/2r) 



(20) 
(21) 



and (pi is given by 

(Pi = iivc sin {di sin d/c) (22) 

In the second case (^1 in 7^ we have 

v'=r'+dl-2rd2 cos (<Po+h) (23) 

where i/^i is given by 

\f/i=SiTC sin (di sin d/c) (24) 

and (/)o is given by (21). 

Once again the formulas for v^ are identical so there 
is no need to make a distinction between the cases. 

Recalling that we are assuming R=l, the situation 
that Q lies inside A is expressed hy u<l or equiva- 
lently u^<l. Similarly ^0 hes inside A if and only 
if ?;^<1. Thus we have a perfectly effective test for 
this situation. 

9.5. Analysis of Cases 

We suppose throughout that both disks B and C 
meet the large circle A, and that B is at least as close 
as C to the center of A. We also assume that the 
radius r of B and C is less than that of A and tliat d 
measures the smaller angle forjned at the center of 
A. These assumptions are equivalent to the alge- 
braic restrictions on parameters contained in (1), 
(2), (3), and (4) of this appendix 1, section 9.2. 

Case Fi'. Both disks are entirely within A. This 
is the case wlien 

c/2 + /'<l. (25) 

This says that disk C is inside A, but we agreed that 
B was at least as close in, so both nmst lie inside. 
We distinguish two subcases. 

Subcase Fn\ B and 6" overlap. This is true when ^^ 

c<2r, (26) 

and the formula for the cojinnon area covered is 

Area=27rr2-|5nC|. (27) 

where \B[]C\ is found according to eq (2). 

Subcase F12: B and C do not overlap. This is true 
when 

c>2r, (28) 



and the formula for the area covered is 
Area = 27^r^ 



(29) 



Case F2: B is entirely inside A, C partially so. 
This is true when rfi + r <l<c/2 + r. Again we dis- 
tinguish two subcases. 

Subcase F21: B and C overlap. The condition is 
(26) and the area is given by 



AYesi=7rr'+\Ar\C\-\BnCl 



(30) 



32 The quantity c denotes the distance of centers of B and C and is calculated 
by c = (rfiHrf22-2rfirf2 cos 0)1/2. 
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where !^4nf| is given by (3), (4), (5) with 6?2 replacing 
c/i in all three formulas for obvious reasons. |^ri C\ is 
given by (2). 

Subcase F22' B and C do not overlap. The condi- 
tion is (28) and the formula for the area is 



Area-TrH+l^nCI. 



(31) 



Case Fs'. Both B and C are only partially inside A 
and they overlap. The condition is (26) and 
rfi+r>l. 

We now consider in order the subcases of F3 
beginning with the case ^— tt and the case c=0. 
We eliminate the possibility c=0 early in the game 
since we call for division by c in many cases. 

Subcase F^q: The centers of B, A, and C are col- 
linear in the order indicated so that the condition is 
^=7r; it can be easily established that in this case 
u<^l necessarily and that the common region BClC 
lies entirely within A. The area is then given by 



Subcase F^i 
tions are 



ATe^=\A[]B\ + \A[]C\-\Bf]C\. (32) 

B and C are coincident. The condi- 
^=0; di=d2 (33) 

and the area is given by 

Area=|^n£|, (34) 

which is calculated by eqs (3), (4), (5). 

Subcase F32 : The centers of A, B, and C are col- 
linear, but B and C are not coincident, and the two 
intersection points of disks B and C lie inside A. 
The conditions are 



^=0; di9^d2; 
where u^ is given by 

u^=r'^-\-did2. 



u'<l 



(35) 
(36) 



It can be readily verified that eq Cl8) reduces to (36) 
under the conditions of (35). The area is given by 



ATeii=^wr'+\Af]C\-\BnC\. 



(37) 



Subcase F^^: The same conditions hold as for 7^32 
except the two intersection points of B and C lie 
outside A. The condition is 



u^=r'+did2>l 
and the area is given by 

Area-i^n5|. 



(38) 



(39) 



In the following situations we shall describe the 
cases according to the conditions on the parameters 
and let the reader figure out the geometry for him- 
self. First we define two new parameters di and 62. 

^,=arc cos (0+d'i-r')l2di) i=l,2. (40) 



Referring to figure 786, ^1 is the angle ^ QA^Bi and 62 
the angle ^CiA^Qo. 

Subcase F^^: The conditions are ^^ 



19^0; u^<l] v''<l', d<d2 
and the area is given by 

Areei=7Tr'+\AnC\-\Br\Cl 
Subcase F35: The conditions are ^^ 

69^0; u'<:i', v'<:i; d>d2 
and the area is given by 

AYe^=\Af]B\ + \Af]C\-\B(]C\. 
Subcase F^^\ The conditions are 
u'^l) v'yi 
and the area is given by 

Ai^Q^=\A{\B\ + \A[\C\-\B[]C[\A\ 



(41) 



(42) 



(43) 



(44) 



(45) 



(46) 



where the last term is the common area of the 
triangular region whose area is calculated in section 
9.3 of appendix 1, eq (14). 

Subcase F^j: The conditions are 



u'<l; 



■■l\ 



d<e2 



and the area is given by 

ATeei=\Af]B\ + \A[]C\-mC[]A\. 
Subcase F^g : The conditions are 

u'<l; v^==l] eyd2 
and the area is given by 

Area-|^n5|+|^nc|-|^nc^|. 

Subcase F^qI The conditions are 

and the area is given by 

Area=:|AnB|. 
Subcase F^q: The conditions are 

and the area is given by 

Area=|An5|. 



(47) 



(48) 



(49) 



150) 



(51) 



(52) 



(53) 



(54) 



33 See section 9.6 of appendix 1 for an explanation of the inequalities on the di. 
See section 9.4 of appendix 1 for definitions of u and v. 

34 0^0 will be true for all cases F^ through Fis. 
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(55) 



(56) 



Subcase F\i: 'Ilie conditions are 

u'' = l; ?;2>1; e>d, 

and the area is given by 

Area=|^n^| + |^nC^|. 
Subcase F42' The conditions are 

u'>v, v'>i; e<e, (57) 

and the area is given by 

Area=:|^n5|. (58) 

Subcase F^^\ The conditions are 

'i/2>i; v^y\] eye, (59) 

and the area is given by 

Kv^^=\A^B\ + \A{\C\, (60) 

The final case occurs when both disks B and C 
meet A but not each other. 
Case Fu' The conditions are 

ch+r>\', rfi+r>l; c>2r (61) 

and the area is given by 

Area^l^n^l + l^nC^i. (62) 

9.6. Special Cases 

Most of the formulas for the area covered that are 
presented in section 9.5 of appendix 1 can be verified 
by a consideration of the geometry of the config- 
urations. An exception to this is the formula 
Trr'^+\A{\C\ — \B^C\ which appears in cases F^2 
and Fu- 

Referring to figure 80 and denoting the areas 
indicated as Gi, G^, and 6^3 we get the following 
expression for the area covered: 



Area-i^nc'i+l^n^i-e^i 

but we also have the equations 

G,^G,= \BW\ 



(63) 

(64) 
(65) 
(66) 



Solving for Gi and substituting in (63) we get 

kx^^=\A[\C\-\B[\C\\'KT\ (67) 

As regards the inequalities on the ^j, if ^<C^2 then 
the portion of disk B that lies outside A is inside C. 
If ^>^2 then this portion does not meet C and the 



area is calculated accordingly. Figures 81a and 
81b refer respectively to cases F34 and F35 and the 
inequalities can be seen geometrically. The inequal- 
ities on e and e, are similarly motivated. 



10. Appendix 



2. Analytical 

iV = 2 



Solution '^' for 



This section deals with the maximization of F{X) 
in the very simple case n=2. As noted in the main 
text, and explained in detail in appendix 1, F{X) is 
given by one of eight different formulas, depending 
on the nature of the configuration formed by the 
fixed circle A of radius i?==l (this is the circle ^Ho be 
covered' ') and the two ^'covering circles'' d and 62 
of radius t<C\. Despite this complication, we shall 
show that the problem can be solved analytically. 

To avoid trivial cases, the assumption l/2<^r<^l 
will be made throughout. As in appendix 1, the 
following notation will be used: 

di=distance from Ci's center to ^'s center, 
c/2== distance from C^2^s center to ^'s center, 
c= distance between 6Vs center and ^2^8 center, 
^== angle between radius of A through d's center 
and that through 62^8 center. 

Thus we have, by the Law of Cosines, 



c''^d,'+d2^-2d42 cos e. 



(1) 



The function to be maximized is given, in set theo- 
retic notation, by 

irF(X)^G{dr,d2,e)=Area (Af\a) 

+Area(^nQ-Area (^nftflft). (2) 

A preliminary remark which greatly simplifies the 
situation is that ^=7r for any configuration which 
maximizes F{X). To prove this, temporarily regard 
di and (4 as fixed, but e as variable. That is, regard 
Ci as fixed but C2 as rotatable around the center ol A, 
Then the first two areas in the right-hand side of eq 
(2) are constant, but the third area is a decreasing 
function of Cand therefore fcf. eq (1)) is a decreasing 
function of for 0< ^< tt and an increasing function 
of e for 7r< e< 2t. 

In what follows, therefore, ^=7r will be assumed, 
so that eqs (1) and (2) become, respectively, 



c=di+d2, 



(3) 



G{dud2,7r)=g{di,d2)=Aresi (AnC,) 

+Area (An Q -Area (ftflQ, ^4) 

where eq (4) follows from the observation that 
AnCir{02=CinC2 when <9=7r. 

Next it will be shown that, for every configuration 
maximizing F(X) , 

di+d2=c<2r, (5) 



di+r>l t=l,2. 



(6) 



35 This solution is due to B. K. Bender and A. J. Goldman (NBS Operations 
Research Section). C. T. Zahn, Jr., suggested several expository improvements. 
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Condition (5) asserts that Ci and C2 meet, while con- 
dition (6) asserts that ft and C2 ^^stick out'^ past A, 
or at least are not entirely interior to A, 

For the proof, first assume tentatively that c>2r 
for some configuration which maximizes F{X). 
Then (U+r<\ cannot hold for both i=l and i=2, 
since this would imply 

4r<c+2r=: {d,+r) + {d2+r)< 2, 

contradicting the hypothesis r>l/2. Thus at least 
one of the lunctions Area {A{]Ci) (i=-l,2) is a strictly 
decreasing function of d^ near the configuration in 
question. For such an i, we can slightly decrease r/, 
and thus increase Area {A^d) without violating 
the condition c>2r. Thus one of the first two areas 
in eq (4) is increased, the other is unchanged, and 
the third remains zero since Ci and C2 are disjoint 
when c>2r. Therefore g{di,d2) has been increased, 
violating the assumption that the original configura- 
tion was maximizing. So the tentative assumption 
that c>2r is untenable, i.e., condition (5) holds for 
every maximizing configuration. 

Now temporarily regard c, and thus the third area 
in eq (4), as fixed, so that eq (3) is a constraint on 
di and ^2. For any configuration in which r/i + r>l 
but d2 + r<il, it would be possible to decrease d^ 
slightly (thus increasing Area {A[\Ci)) and to in- 
crease ^2 by the same amount so that d2+r<l is not 
violated and Area {A[]C2) retains the value Trr\ 
Thus g(di,d2) would be increased, and so the original 
configuration could not have been maximizing. A 
similar argument apphes with i=l and i=2 inter- 
changed. Therefore a maximizing configuration 
either obeys f6) for i=l,2, or obeys 



di+r<l iori=l,2. 



(6a) 



Under the condition (6a), however, the first two 
areas in eq (4) have the value xr^ and only the third 
one IS variable. This area is minimized (i.e., 
gidiA) is maximized) by choosing d^ and ^2 (and thus 
c=-di+d2) as large as possible. Subject to (6a) 
these choices are di-l-r(i=l,2), which still satisfy 
(5) since 

c=2—2r<2r because r>l/2. 

But these choices also obey (6). This completes the 
proof that (5) and (6) hold for all configurations 
maximizing F(X). 

In what follows, therefore, conditions (5) and (6) 
will be assumed. It is convenient to introduce the 
following quantities : 

2<9i=- angle intercepted at ^'s center by subtended 

arc of d, 
2 (^, = angle intercepted at C/s center by subtended 
arc of A, 
^i=length of common chord of A and Cf, 
2\[/=eing[e intercepted at center of either d or ft 
by subtended arc of the other of ft or ft, 
2=- length of common chord of ft and ft. 



It is readily found that 

cosd,=:{l-r'+d,')/2d,, (7) 

cos6i=(l-r'-d,^)/2rd„ (S) 

cos \p=c!2r (9) 

sin ei = r sin 0^ = ^3,., r sin \p=\z, (10) 

Area (Anft-) = 7rr2+6/,— i sin 2(9,-rVz + ir2 sin 2(^,, 

(11) 
Area (ftnft)=2rV-r' sin 2;^. (12) 

From the geometr}^ of the situation, it follows that 
b (Area {A[]Cd)IHdi) = -Zu (13) 

an analytical derivation of this will be given later. 
Exactly the same argument shows that 

5(Area (ftnft))/dd,— -2, 

so that (see eq (4)) we have 

<)glM, = z-z, (^=l,2). (14) 

Now it will be shown that there is precisely one 
maximizing configuration, the one characterized 
rather elegantly by 



Zi=22=^Z 

or equivalently, via eq (10), by 



(15) 



(16) 



Note that, as might be expected, the maximizing 
configuration is symmetric in the sense that di=d2. 

To prove eq (15), tentatively suppose it false. 
Without loss of generahty suppose Zi^^z. Then the 
function g{di,d2) assumes its maximum on the tri- 
angle : 



T: rfi>l-r, d2>l- 



d^ + d2<2r, 



defined in the (di,d2)— plane by conditions (5) and 
(6), at a point at which dg/ddi does not vanish. 
Such a point must lie on the boundary of T, and in 
fact not on its horizontal leg (endpoints excluded) ; 
this follows from standard calculus arguments. If 
the point is on the vertical leg of T (upper end-point 
excluded), then on the one hand dg/ddi must be non- 
positive and thus negative, so that z<Czi by eq (14), 
and on the other hand di+r= 1 so that ft is internally 
tangent to A, implying Zi=0. Since z<Czi and 
01=0 are incompatible, this case is ruled out. If 
the point is on the hypotenuse of T ("upper endpoint 
excluded), then on the one hand dg/ddi must be non- 
negative and thus positive, so that 2>0i, and on the 
other hand c=2r so that ft and ft are internally 
tangent, implying z=0. Since zy>Zi and z=0 are 
incompatible, this case is also ruled out. Finally, if 
the point is the upper vertex of T, then on the one 
hand consideration of the directional derivative 
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along the hypotenuse of T yields 



-^i<0, (17) 



while on the other Jiaiid ^i = and ^ = as above. 
This imphes ^2=0, so that Ci and G are externally 
tangent to each other and internally tangent to A. 
Such a configuration can only occur if r=l/2, con- 
tradicting our assumption that ?'>l/2. So every 
alternative to eq (15) has been ruled out, and the 
equation must liold. To describe the maximizing 
configuration more explicitly, let x denote the com- 
mon value of di and ck. Equations (8), (9), and (16) 
then yield 

which implies that 

x-((l-r2)/3)^/2 (jg^) 

As a check, note that x— >l/2 as r-^1/2 and x-^0 as 
r-^1, as would be anticipated. 

The maximuni coverage ratio T^max can now be 
found in terms of r. First we have, from eqs (4), 
(11), (12), and (16) 

7rF^ax = 27rr2+(2^,-4rV)-sin 2^,; + 2r2 sin 2^. 
From eqs (7) and (18), however, 

cos ei=2{{l-r')l3y^'=2x, 
which with the aid of eqs (9), (10), and (16) yields 
2r^ sin 2\l/=4r^ sin <pi cos \l/=4:r(su\ 6i){x/r) = 
2 sin di cos ^i=sin 26 i. 

Thus the last two terms in the above expression for 
TrFmax cancel each other, leading to 

7rF„iax=27rr2+2arc cos (2a:)-4r2arc cos (x/r), (19) 
where x is given by (18). 



Table 11. Comparison of analytic and computer solutions'" 
for the case n = 2 



r 


X 


Xc 


Fmarlir 


RATI0„,„.. 


?46 - - - 


0. 47735 
. 45069 
.41926 
.38188 
. 33657 
. 27951 
. 20091 


fc 0. 47754 
. 450()9 
.41753 
.38181 
. 33802 
. 28128 


0. 601 
. 68B 
. 762 
.829 
.889 
.939 
.979 


0.600 


5^ 


. 68i) 


iHe 

3/4 

me 

^^ 

l>)i6 


.762 
.829 
.889 
.939 
.979 







« X is the distance of the centers of each of C\ and C2 from that of . 1 in t he opt inial 
configuration, and Fmaxjir is the value of coverage obtained from t he configura- 
tion, i.e., the maximuni coverage. The corre.sponding values .tc and liiWiOmai 
are those obtained from the computer simulation at a mesh of 25(). 

^ At the coarser mesh of 64, the value is Xc=^A'l&. 



The value of x and maximum coverage (i^max/^r) 
are compared in table 11 with the corresponding 
values Xc and EATIOmax obtained by the computer .^^ 
As can be seen, the agreement is excellent. 

We conclude with an analytical derivation of ec^ 
(13). 

First use eq (11) to wtite 

d(Area {A[]Cd)IWi)='^ sin^ ^,(d^,/d((/,)) 

-2r2sin2^,(5^,/d(rfO). 
By eq (10), tliis can be written 

d(Area {A[]Ci))lb{dd = Zi shi ^,(d^,/d(r/,)) 

-rsin^,(d^,/5(r/0). (20) 

From eqs (7) and (8), however, 

-sin ^,(d^,/d(./,)) = - {l-7''-dt')l2d,\ 
-r sin ipi{biPilb{dd)=^-il-r'+di')l2di\ 

Substitution of these results into eq (20) vields 
eq (13). 



39 This table was prepared by C. T. Zahn, Jr. 
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11. Appendix 3. The FORTRAN Program 



C CYCLING FIRST POSITIVF GRADIENT 

-S-nKE-OCT— 1 onnrLQ_a 

S P2 OCT 237105624673 
RFAD1001 



PRiNTinoi 
1001 FORMAT( 72H 



1 ) 

DIMEN 51 ON A(50 )» IP(30) »JQ(50)»IP0(50) f JQ0(50)»K2(1000)»IT( 4 »51:?)' 

XJT{4,512)fIB{4«512)fJRl4»512) 

DIMENSION NX(4) »IPX(50) »JQX(50) 

RFADltN7 

IDX2 = 1 ^ 

CALL MK2T(K2) 

PRINT 2010 

1000 READlfJ7»!<0»NC»N»KBlGtJXPOf JFACT»IR0»N4 

1 F0RMAT(7I10) 



IDX=1 
CLA P2 



STO 14 

PRINT 2011»IDX2 



PRINT 2012fNC»IR0»K0iN»KBIGiN4 
K=:K0 



IR^IRO 
N F W X=0 



JXP=JXPO 
50 IF( J7)3f lOtlO 



10 READ9* (IP (IX)* JQ(IX)» IX=lfNC) 
GO TO 4 



3 CALL XBAR( JXPfKf IP*JQfK2#NCf 14) 

4 CALL SUMX(K»IP fJQ i I R » NC »N SUM » K2 ) 



PRINT 2020f IDX*( IX»IP( IX)»JQ( IX) f IX = 1»NC) 
CALL RAT (K»K2fNSUM»RATI0) 



PRINT 2021* RATIO 
PRINT 2022» K 



KS0=K2(K) 
IRQ=K2( IR) 



CALL VECTOR( iTf JTflBf JB»K2»IR) 
CLA IR 



S ALS 1 
S STO IR3 



S SUB ONE 
S STO MAX 



9 FORMAT( 1017) 

7 IF( I0-NC)101f 102*102 



102 IX4=NC 
10 = 1 



GO TO 201 
101 1X4=10 



10=10+1 
GO TO 201 



201 CALL NEAR ( NS • I * I PX ♦ JQX » I R 3 t IP t JO iNC ) 
INT5 = IP( 10) 



INT6=J0( 10) 
INT7=K-IR-1 



IF(K2(INT5)+K2( INT6)-K2( INT7) ) 1234,1234»1235 

1234 IF(NS) 1238»2005*1238 

1238 DO 1300 IXS=1*4 
CALL SCANl (IT*IB*JT»JB» I0»K2»JD»IP»JQ»IXS»MAX»NSf lPX»JQXfKSQ»IRQ) 

IF( JD) 1300»1300*1307 
1300 CONTINUE 

GO TO 2005 

214 



1235 DO 1?^0 lXS = li4_. 



CALL SCAN? ( TT» TBOTfJB* I0,K2 ♦ JD ♦ I P ♦ JO , I XS ♦ MAXi NS» IPX t JQX i KSQ • I RO) 
IF( JD) l?^0»125O»1307 , 



1250 CONTINUF 

2 005 IF( I0~IX4) 304»?n^»302 



30? 
305 


IF( 10- 
10 = 1' 


-NC)304! 


i305 


^305 


304 


GO 
10 = 


TO 
:I0- 


201 
fl 






303 


GO TO 
IF(K- 


201 
CBir.) IR 


n^t 


19 



1307 1X0=10 
MX2=IXS 



CALL XNFW (MX2»IPf JQf IXO) 

NFWX = NEWX-H 

GO TO 7 
18 PRINT 2027> NFWX 



PRINT 2030tKf(IX5lP(lX)fJQ(IX)fIX*l»NC) 
CALL 5UMX ( K, IP»J0«IRiNCfN.SUM»K2 ) 



CALL RAT(KfK29N5UMfRATIO) 
PRINT 2 021 f RATIO 

"CAIL 'refine' ("ITVjT* IR»JB«NCiK«IR#IPiJQf JXPf JFACT)' 

PRINT 2032» K 

GO TO 6 
_1_9 CAL_L 5UMX (K»JP_*JQ»IRtNC»NSUM»K2 ) 

CALL"RAT(KfK2\NSUM7RATI0) 
P RINT 2040» K »( IX»IP( IX) ♦ JQ ( IX) » JX«1»NC) 

PRINT 2041f RATIO 

PRINT ^04?y NFWX 

53 IF( IDX-N4)41 »42»42 

41 IDX=IDX^-1 



GO TO 2 
4? IF( IDX2_-N7)43jl4Ai44 
44 CALL RETURN 
43 IDX2=IDX?-H 



GO TO lOOO 

2010 FORMAT( IHl » l_?Xj_38HCYCL I_NG _F_I_B3T POSIT I VE GRADIENT SEARCH) 

2011 FORMAT{lHlfl9Xf5HCASE 112//) 

2 012 FORMAT(2OXfl0HTHFRF ARE •I2f26H COVERING DISCS OF RADJ US^ f H* IH/ ♦ I 
X3/20X»25HMONTE CARLO IS USED WITH #I3t8H T R I AL S . /20X f 19HF INAL mESH 
X SIZE IS »I4»13H, KICK OFF »I3»7H TIMES^) 

2020 FORMAT! ///15X*22HINITIAL CONFIGURATION , I 4 // ( 1 5X » I 2 # 4H ) ♦I5fl8)) 

2021 F_PRMAT ( // 15 X »QHR AT I IS fFlO.e) 

2022 F0RMAT(/////10X*8HMESH IS 114//) 

2027 FORMA T! //10X»I3>21H MOVES HAVE BEEN MADE) _ 

2030 F0RMAT(7//i5X9 3lHRELATIVE MAXIMUM UNDER MFSH OF fT4»21H'rS THE CON 

1FIGURATI0N//(15X5I2»4H) ♦ I5»I8 ) ) 

2032 F0RMAT(/////10X»19HMESH IS REFINED TO fI4) 

2040 F0RMAT(7///15X#41HRELATIVE MAXIMUM WITH FINAL MESH SlZE_0F fI4,21H 
X IS, THE C0NFIGURATI0N//(15X#I2,4H) #15,18)) 

2041 FOR MAT! //15X»2 4HFINA L VAL UE OF RATIO IS_iFj0._6J 

2042 FORMAT(7/15X»26HTOTAL NUMBER OF MOVES WAS ,1 3/////5 (lOX • lOHXXXV XXX 
XXXX )///// ) ' 

END 
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